set more off


local ID `1'
local vij `2'
local label `3'
local J `4'
local altID `5'

*Use more robust logit probability calculation
gen double sij`label' = .
forvalues j = 1/`J' {
gen onlyj = `vij' if altID == `j'
bys ID: egen maxvij = max(onlyj)
gen expvik = exp(`vij'-maxvij)
bys `ID': egen sumexpvik = sum(expvik)
replace sij`label' = 1/sumexpvik if altID == `j'
drop onlyj maxvij expvik sumexpvik
}


